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Some of the most interesting scenarios that can be studied in astrophysics, contain fluids and 
plasma moving under the influence of strong gravitational fields. To study these problems it is 
required to implement numerical algorithms robust enough to deal with the equations describing 
such scenarios, which usually involve hydrodynamical shocks. It is traditional that the first problem 
a student willing to develop research in this area is to numerically solve the one dimensional Riemann 
problem, both Newtonian and relativistic. Even a more basic requirement is the construction of the 
exact solution to this problem in order to verify that the numerical implementations are correct. 
We describe in this paper the construction of the exact solution and a detailed procedure of its 
implementation. 



PACS numbers: 04.40.-b,04.25.D-,95.35.+d,95.36. 
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I. INTRODUCTION 

High energy astrophysics has become one of the most 
important subjects in astrophysics because it involves 
phenomena associated to high energy radiation, modeled 
with sources travchng at high speeds or sources under the 
influence of strong gravitational fields hkc those due to 
black holes or compact stars. Current models involve a 
hydrodynamical description of the luminous source, and 
therefore hydrodynamical equations have to be solved. 

In this scenario, due to the complexity of the system of 
equations it is required to apply numerical methods able 
to control the physical discontinuities arising during the 
evolution of initial configurations, for example the evolu- 
tion of the front shock in a supernova explosion, the front 
shock of a jet propagating in space, the edges of an accre- 
tion disk, or any shock formed during a violent process. 
The study of these systems involve the implementation of 
advanced numerical methods, being two of the most effi- 
cient and robust ones the high resolution shock capturing 
methods and smooth particle hydrodynamics which are 
representative of Eulcrian and Lagrangian descriptions of 
hydrodynamics, each one with pros and cons. 

It is traditional that a first step to evaluate how ap- 
propriate the implementation of a numerical method is, 
requires the comparison of numerical results with an ex- 
act solution in a simple situation. The simplest problem 
in hydrodynamics is the ID Riemann problem. This is 
an excellent test case because it has an exact solution in 
the Newtonian case (e.g. [l|) and also in the relativistic 
regime [1, |1] , where codes dealing with high Lorentz fac- 
tors are expected to work properly. From our experience 
we have found that the existent literature about the con- 
struction of the exact solution is not as explicit as it may 
be expected by students having their first contact with 
this subject. This is the reason why we present a paper 
that is very detailed in the construction and implemen- 
tation of the solution. We focus on the solution of the 



problem and omit some of the mathematical background 
that is actually very well described in the literature. 

The paper is organized as follows. In section [11] we 
present the Newtonian Riemann problem and how to im- 
plement it; in section Hill we present the exact solution to 
the relativistic case and how to implement it. Finally in 
section ITVl we present some final comments. 



II. RIEMANN PROBLEM FOR THE 
NEWTONIAN EULER EQUATIONS 

The Riemann problem is an initial value problem for 
a gas with discontinuous initial data, whose evolution is 
ruled by Euler's equations. The set of Euler's equations 
determine the evolution of the density of gas, its velocity 
field and either its pressure or total energy. A comfort- 
able way of writing such equations involves a flux balance 
form as follows 



dtxv + a,F(u) = 



(1) 



where u = (mi, ti2, us)"^ = (p, pu, E)^ is a set of conserva- 
tive variables and F is a fiux vector, where p is the mass 
density of the gas, v its velocity and E = p{^v^+e), with 
£ the specific internal energy of the gas. The enthalpy of 



the system is given by the expression H = 



-h, where 



h is the specific internal enthalpy given by h — e + p/p^ 
where p is the pressure of the gas. The fiuxes are ex- 
plicitly in terms of the primitive variables p,v,p and the 
conservative variables [l[ 



F(u) 
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5(3 -r)g + (r- 1)1.3 



The initial data of the Riemann problem is defined as 
follows 
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\ Uu, X > Xq, 

where and uj^ represent the values of the gas prop- 
erties on a chamber at the left and at the right from an 
interface between the two states at x = xq that exists 
only at initial time. 

The evolution of the initial data is described by the 
characteristic information of the system of equations, and 
this is why the properties of the Jacobian matrix are im- 
portant. The Jacobian matrix of the system of equations 
is ^(u) — ^ and explicitly reads 



A- ^(r-ay (3-r)t- r-i . 

Its eigenvalues satisfy the condition Ai(u) < A2(u) < 
A3(u) and are given by 



combination {x — xo)/t] it can be seen that such behavior 
implies that the following conditions hold Q 

dui du2 du3 
r\ r| r| 

where i indicates the component of a given eigenvector. 
On the other hand, the Rankine Hugoniot conditions re- 
late states on both sides of a shock wave or a contact 
discontinuity 

AF = VAu, (6) 

which arc simply jump conditions, where Au is the size 
of the discontinuity in the variables, V is the velocity of 
either the contact discontinuity or shock and AF is the 
change of the flux across the discontinuity. 

A. Contact discontinuity waves 



a 



Xi = V 
X2 = V 
A3 = v + a 



(2) 
(3) 
(4) 



The contact discontinuity is described by the second 
eigenvector and evolves with velocity A2. Let us then 
analyze the second eigenvector. In this case the Riemann 
invariant conditions read 



where 



is the speed of sound in the gas, which 



depends on the equation of state. For the ideal gas p = 
p£(r — 1), where F is the ratio between the specific heats 
at constant pressure and volume F = Cp/cv, the speed of 

sound is a = ^J'^■ On the other hand, the eigenvectors 
of the Jacobian matrix read 





The eigenvectors ri , Y2, are classified in the following 
way: 

• they are called genuinely non-linear when satisfy 
the condition V„Ai • ri(u) 7^ 0. 

• and linearly degenerate when V„Ai • Vi{u.) = 0. 

It happens that r2 is linearly degenerate and repre- 
sents a contact discontinuity, however the other two are 
genuinely non-linear. 

Depending on the particular region of the solution we 
will use both the Riemann invariant conditions for rar- 
efaction waves and the Rankine Hugoniot conditions for 
shocks and contact discontinuities. The Riemann invari- 
ants are based on the self-similarity property of the so- 
lution in some regions, in the sense that the solution de- 
pends on the spatial and time coordinates (z, t) with the 



dp d{pv) dE 
1 V 

These relations implies that d{pe) = dv = 0, further 
implying that p = constant and v = constant across 
the contact wave. In order to relate the two sides from 
the contact discontinuity we use the Rankinc-Hugoniot 
conditions, which are given by 



Plvl - Prvr = 
PlvI+pI - PRvji+p% = 
vlIEl +Pl) - vrIEr+pr) = 



VcipL-pR), (7) 
VcipLVL - Prvr),{^) 
Vc{vl{El+Pl) 
VRiER+PR)). (9) 



Here Vc is the velocity of propagation of the contact dis- 
continuity. 

The discontinuity travels at speed X'^ — v therefore 
the Vc = V. For this reason from equation (O follows 
that vl — vn = Vc- As a consequence of this, equation 
dHl) gives the condition pL = Pr, which implies ^ is 
satisfied. Notice that no condition on the density arises, 
which allows the density to be discontinuous. 



B. Rarefaction waves 

At this point we do not know the nature of waves 
1 and 3, and we can assume they may be rarefaction 
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waves. Once again we use the Riemann invariant equali- 
ties, which for vectors 1 and 3 read 



dp d(pv) dE 

1 V — a H — av^ 

dp d{pv) dE 

1 V + a H + av 

Manipulation of these equalities results in the following 
equations 



dp 

dv 
dp 

dv 
de 
dp 



P 
a 
V_ 
p2 



forAi, (10) 
for A3, (11) 
for both Ai and A3. (12) 



The next step is to integrate these equations assuming 
an equation of state, in our case the ideal gas. From p2)) 
we obtain 



■p = Kp 



(13) 



where K \b a, constant. A rarefaction process is isen- 
tropic (unlike a shock), and therefore the states at the 
left and at the right from the wave obey with the 
same constant K. 

Using this expression for p in the speed of sound we 
have a = \J KV (F^^ = ^JV^p|p^ which substituted into 
(jlOllip rcsuhs in 



v = ± J VKTp^-^dp + k = ±^^^ + k, 



(14) 



where + stands for the wave moving to the right (the 
case of A3 and corresponding to a rarefaction wave) 
and — when moving to the left (the case of Ai and ri 
corresponding to a rarefaction wave), where k is an in- 
tegration constant and therefore the velocity is constant 
as well. This property allows us to set relations between 
the velocity of the gas on the state at the left and at the 
right from the rarefaction wave, explicitly there are two 
possible cases: 

i) When the wave is moving to the left, condition 
implies that 



2ai_ 

r- 1 



= VR 



r- 1' 



(15) 



ii) When the wave is moving to the right, condition 
(fT4|l implies 



2ai_ 

r- 1 



2aR 

r - 1 



(16) 



When the wave is moving to the left, we assume in- 
formation from the left state is available and we look for 
expression of the variables on the state to the right from 
the wave. For the velocity of the fluid at the right state 
we then have from (fTSl) 



VR = VL 



r- 1 



[an. - cll] 



(17) 



now considering that the speed of sound on both sides 
obeys a = y^KTp^ = ^/Tp/p (see ([T^) 



PR 

PL 



(18) 



a useful expression for vr arises 



VR = VL 



2aL_ 

r- 1 



PR 

PL 



(19) 



The only unknown quantity is pR. 

On the other hand, when the wave is moving to the 
right we assume we know the information at the state at 
the right from the wave, then we search for expressions 
of the variables on the state at the left. For the velocity 
we find according to 



VL = Vr 



r- 1 



[aR - o-l] 



(20) 



and the speed of sound on both sides obeys 



a-L = aR 



Pl_ 
PR 



(21) 



which finally implies 



VL = Vr 



'2aR 

r- 1 



1 - 



Pl_ 

PR 



(22) 



The only unknown quantity in this case is pl- 

The rarefaction zone has a finite size, bounded by two 
curves, the tail and the head. The head of the wave is the 
line of the front of the wave and the tail is the boundary 
left behind the wave. The region in the middle is called 
the fan of the rarefaction wave. 

The velocity of all the particles between the head and 
the tail obeys the following expression 



Xq 



t 



V ± a, 



(23) 



where -I- is used when the wave is propagating to the right 
and the — when it is moving to the left. Then, when the 
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wave is moving to the left, using this expression we have 
Oii = 'VR — {x—Xo)/t, which substituted into ([T5)) provides 
the following expression for the velocity of the gas on the 
state at the right from the wave is 



1/-P .^ , x-xq 



(24) 



Then it is possible to calculate the pressure and density 
as well. Substituting (|24p into ((TS|) and ((T5)) we obtain an 
expression for the pressure also at the state to the right 



PR = PL 



r - 1 



r + l aL{T + l) 



VL 



X — Xq 
t 



(25) 

Now, using this into ([T3| implies the expression for the 
density 



C. Shock waves 

Similar to the previous case, the shock can move either 
to the right (if A3 and correspond to a shock wave) or 
to the left (if Ai and ri correspond to a shock wave) , and 
for each of the two cases there is known and unknown 
information. When a shock is moving to the right one 
is expected to have information of the state at the right 
from the shock and conversely, when the shock is moving 
to the left one accounts with information of the state at 
the left. 

Shocks require the use of Rankine Hugoniot conditions 
(|ni). We express these conditions in terms of the primitive 
variables as follows 



PlVl - Prvr 
PR 

vl{El +Pl) - vr{Er +Pr) 



PLvl +PL- PRVr 



S{pL - Pr), 
S{plvl - Prvr), 
S{El-Er), 



PR = PL 



r- 1 



r + i aL(r + i) 



VL 



X — Xo 

t 



(26) 

Then finally we have expressions for the velocity, pressure 
and density on the state at the right when the wave is 
moving to the left. 

Similarly when the wave is moving to the right we have 
from ((23|) that ul = vl + (x — xo)/t, which substituted 
into (|22p implies the following for the velocity on the 
state at the left from the wave 



VL 



r + l 



, ^f-r 1^ x-xo 
-o-R +2' '^''^ 1 — 



(27) 



In order to obtain the expressions for the pressure and 
the density, we substitute this last expressions into (|16p 
in order to relate the speeds of sound, and then using 
(|2ip we finally obtain the expression for the pressure at 
the left 



PL =Pr 



r- 1 



r + l afl(r + i) 



VR 



X — Xo 
t 



(28) 



Finally using the equation ()13|) we obtain the density 



where S is the speed of the wave, which may take the 
values v—a or v+a depending on whether the wave moves 
to the left or to the right respectively. Manipulating these 
equations one gets 



Plvl 



Prvr, 



PlvI+pl = PRvj^+PR, 
Vl{El+Pl) = Vr{Er+Pr) 



(30) 
(31) 
(32) 



where vl = vl — S, ~ vr ~ S are velocities in the 
rest frame of the shock and El — PL {\v\ +£l) and 
Er = PR (^Wfl + £r)- These expressions correspond to 
the Rankine Hugoniot jump conditions measured by an 
observer located in the rest frame of the shock wave. 

From equation ([5(1)) , we introduce the mass flux defini- 
tion 



J = plvl = PRVR. 



(33) 



Then, from equation (|3T|) and the mass flux definition 
before mentioned , we can get an expression for j, which 
is given by 



J 



PR -PL 

VR - VL 



PR- PL 
vr-vl' 



(34) 



PL = PR 



r- 1 



r + l aR{r + i) 



VR 



X — Xq 
t 



(29) 

In this way we have relations between the variables on to 
the state at the left and at the right from a rarefaction 
wave. These relations will be useful when solving the 
Riemann problem. 



which is a consequence of j being invariant under 
Galilean transformations. Considering the shock is mov- 
ing to the left, we would be interested in constructing the 
variables on the state at the right from the shock and we 
can start with the velocity, which can be written as 



VR = VL 



PR -PL 
j 



(35) 
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Now, in order to express the velocity in terms of the 
pressure and the variables of the state at the left from 
the shock, we can rewrite (l33l) as follows 



(36) 



PR PL 

Thus, substituting this into p4p we obtain 
• 2 PR - PL 



1 

PR 



1 

PL 



(37) 



On the other hand, using equation ([32]) and the expres- 
sion for the specific internal enthalpy h we can easily 
get the following expression for the difference of internal 
specific enthalpies 



1 



(38) 



where /il = £l +Pl/pl and hn = er + pr/ pr. Now, 
from equations (|30p and pip we give expressions for the 
vclocititcs measured by the observer located in the rest 
frame of the shock wave 



PL PL - PR 



PR PL 
PR PL 



PR 
PR 



PL PL - PR 



With the substitution of these last equations into ((38|) 
and considering the definitions for the specific internal 
enthalpy mentioned above, we obtain 



£r - El - - 



1 (PL +PR.){PR - Pl) 



PL PR 



Assuming the gas obeys an ideal equation of state we get 
an expression for the density as follows 



PR ^ PL(r-i) + Pfl(r + i) 
PL PR{T-i)+pL{T + \y 



(39) 



Notice that this expression relates the density among the 
two sides from the shock. Now, substituting this expres- 
sion into (1371) we obtain 



r - 1 
Bl = TT—rPL- 



f = '-^. Al 



Al 



(r + i)pL^ 



r + i 



(40) 

Thus, the expression for the velocity psp can be written 
as follows 



VR = VL - (pr-Pl)] 



A, 



Pr + Bl 



(41) 



From expression ((36|) and using (|40| we express the shock 
velocity as follows 



S = I'L 



/pfl(r + i)+pi(r-i) 



2pL 



Finally, using the sound speed expression = \J^^ we 
obtain the final expression for the shock velocity 



VL - aL\ 



^ [T + 1)PR 

2plT 



2T 



(42) 



Analogously, when the shock moves to the right, it is 
possible to construct the expressions for the variables for 
the state at the left from the shock 



VL = VR + {p 



PR) 



Ar 



PL = PR 



^ PL+Br' 

PRjT - 1) + pLjT + 1) 
Pl(T-1)+pr{T + 1)' 



vr + aRi 



'(r + i)pL 



1 



(43) 
(44) 

(45) 



, ^prT 2T 
and we let this as an exercise to the reader. 

D. Classical Riemann Problem 

The Riemann problem is physically a tube filled with 
gas which is divided into two chambers separated by a 
removable membrane at x = xq. At the initial time the 
membrane is removed and the gas begins to flow. Once 
the membrane is removed, the discontinuity decays into 
two elementary, non-linear waves that move in opposite 
directions. 

Depending on the values of the thermodynamical vari- 
ables in each chamber, four cases can occur. Considering 
the fluid is described on a one-dimensional spatial do- 
main, rarefaction and shock waves can evolve toward the 
left or right from the location of the membrane. 

In general the solution in all the cases can be studied 
in six following regions: 

Region 1: initial left state that has not been yet 
influenced by rarefaction or shock waves 

Region 2: wave traveling to the left (may be rar- 
efaction or shock) 

Region 3: region between the wave moving to the 
left and the contact discontinuity, called region 
star-left 

Contact discontinuity 
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Region 4: region between the contact discontinuity 
and the wave moving to the right, called region 
star-right 

Region 5: wave traveling to the right (may be rar- 
efaction or shock) 

Region 6: initial right state that has not been yet 
influenced by rarefaction or shock waves 

Regions 2 and 5 are special. If the wave propagating in 
such regions is a rarefaction wave the region involves a 
head-fan-tail structure, whereas if it is a shock the region 
becomes only a discontinuity. Counting from left to right 
on the spatial domain, the results can be reduced to the 
following four possible combinations of waves: 

1) rarefaction-shock 

2) shock-rarefaction 

3) rarefaction-rarefaction 

4) shock-shock 

with a contact discontinuity between the two waves in 
all cases. It is worth noticing that these combinations 
can occur under a wide variety of possible combinations 
of the initial values of the thermodynamical variables. 
In this paper we illustrate each of these scenarios using 
particular sets of initial conditions. 



Region 3 plays the role of the state at the right from 
the rarefaction wave and region 1 the state at the left. 
Then we can use ()19|) to obtain an expression for 



V3 = Vi 



2ai 

r- 1 



^1 -1 
pi 



(46) 



On the other hand, region 4 plays the role of a state at 
the left from the shock wave and region 6 the role of the 
state at the right. Then we use (|43p to calculate v^: 



1)4 



V6 + (P4 -P6)\ 



A. 



Vi + Be 



(47) 



where Ag = 2/pe/{T + 1) and Bq = peiT - l)/(T + 1). 
Given that V3 — — v* , equating both expressions one 
obtains a trascendcntal equation for p* : 



A. 



Be 



2ai 

r- 1 



+ve — vi = 0. 

(48) 

Unfortunately as far as we can tell, no exact solution is 
known for p* , and then we proced to construct its solution 
numerically. Once this equation is solved, ps and pi are 
automatically known, and ^3 and 1)4 can be calculated 
using (|46|) and (|47)) respectively. 

Then, it is possible to calculate ps using at both 
sides of the rarefaction zone, given C is the same on both 
sides because it is an isentropic process: 



1. Case 1: Rarefaction-Shock 

This case corresponds to the typical case used to test 
numerical codes, a test called the Sod's shock tube prob- 
lem . A traditional set of initial values that produces 
this scenario corresponds to a gas with higher density and 
pressure in the left chamber than in the right chamber, 
and the velocity is set initially to zero in both. 

A rarefaction wave travels into the high density region 
(moves to the left), whereas a shock moves into the low 
density region (moves to the right). 

Summarizing, the problem then involves five regions 
only. Regions 1 correspond to the initial state to the 
left that has not been influenced by the evolution of the 
system. Region 2 corresponds to a rarefaction wave con- 
taining the head-fan-tail structure, region 3 and 4 are the 
left and right states separated by the contact discontinu- 
ity. Region 5 reduces to the shock. Finally region 6 is 
the initial state at the right chamber that has not been 
influenced by the evolution of the system. 

The goal is to determine the state in all the regions 
using the relations between the thermodynamical quan- 
tities constructed before. 

The starting point to construct the solution happens at 
the contact discontinuity, where the velocity and pressure 
obey the conditions ps = Pi = p* and V3 = V4 = v* . 



pa 



, P3 
Pi — 



i/r 



(49) 



where now pi , pi and p^ are known. On the other hand 
one can also calculate p4 using 



P4 = Pe 



P6(r-i) + P4(r + i) 
P4(r-i) + p6(r + i) 



(50) 



also in terms of known information. With this informa- 
tion it is already possible to construct the solution in the 
whole domain. We explain how to do it region by region. 
A scheme of how the regions are distributed is shown in 
Fig. [II 

1. Region 1 is defined by the condition x—xq < tVhead, 
where Vhead is the velocity of the head of the rar- 
efaction wave given by the characteristic value of 
the Jacobian matrix evaluated at the location next 
to the head from the left side, that is, considering 
([2]) = ui — fli. The solution there is simply 



Pexact — 7^1 1 
^exact — ^^1 1 
Pexact — Pi- 
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Region 2 is defined by the condition tVhead < 
X — Xq < tVtaii , where Vtau is the same characteris- 
tic value again, but this time evaluated at the tail 
curve, that is VtaU = V3 — a^. This is the fan region 
for a rarefaction wave moving to the left, for which 
we simply use expressions ()24|25|26|) that need only 
information from region 1 and obtain 



Pexact 



Pexact 



'^exact 



Pi 



Pi 



2 r- 1 
r + i ^ ai(r + i) 

2 r - 1 



r + i ai(r + i) 



t 



X — Xa 



t 



r + 1 



1 s X ~ Xr) 



(51) 



3. Region 3 is defined by the condition tVtaii < x — 
Xq < tVcontact, where Vcontact is the velocity of the 
contact discontinuity, which is the second eigen- 
value ([3]) of the Jacobian matrix evaluated at this 
region, that is Vcontact = V3 = V4. The solution 
there finally reads 



Pexact — PS 1 
'^exact — ^3; 
Pexact — PZ- 



Case 


PL 


VR 


VL 


VR 


PL 


PR 


Rarefaction- S ho ck 


1.0 


0.1 


0.0 


0.0 


1.0 


0.125 


Shock- Rarefaction 


0.1 


1.0 


0.0 


0.0 


0.125 


1.0 


Rarefaction-Rarefaction 


0.4 


0.4 


-1.0 


1.0 


1.0 


1.0 


Shock-Shock 


0.4 


0.4 


1.0 


-1.0 


1.0 


1.0 



TABLE I: Table with the initial data for the four different 
cases. We choose the spatial domain to be a; £ [0, 1] and the 
location of the membrane at xo = 0.5. In all cases we use 
r = 1.4. 



1 










6 












1 






2 




3 


4 


6 



FIG. 1: Description of the relevant regions for the 
Rarefaction-Shock case. 



4. Region 4 is defined by the condition tVcontact < x — 
Xq < tVshock, where according to (|45|) . the velocity 
of a shock moving to the right separating regions 

4 and 6 is Vshock = ve + «6 ^^2rp!f + where 
06 = \/P6^/ Pe- Then the solution in this region is 

Pexact = 7*4, 
Vexact = Vi, 
pexact = Pi ■ 

as calculated 

5. There is no region 5. 

6. Region 6 is defined by tVshock < x — Xq. In this 
region the solution is simply 



2. Case 2: Shock-Rarefaction 

This case is identical to the previous one, except that 
we choose the initial pressure and density are higher on 
the right chamber. After initial time, the wave traveling 
to the left is a shock, while the one moving to the right is 
a rarefaction wave. This implies that region 2 plays the 
role of region 5 in the previous case and region 5 has the 
tail-fan-head structure of a rarefaction wave. 

Starting from the contact discontinuity, the conditions 
V3 ~ V4 ~ V* and Ps ~ Pi ~ p* hold. The conditions on 
a shock wave moving to the left imply according to (|4ip 
that the velocity of the state at the right is 

y,=y,-ip,-p,)^l^l, (52) 



Pexact P6i 
exact ^6: 
Pexact — PQ • 



An example of how the solution looks like is shown in 
Fig. Hfor initial data in Table H] 



and information from the rarefaction wave interface can 
be obtained from (P7)) for V4 as the velocity on the state 
at the left from a rarefaction wave moving to the right 



V4 = Vq 



2ae 



(53) 
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FIG. 2: Exact solution for the Rarefaction-Shock case at time 
t — 0.25 for the parameters in Table H] 



Equating these two expression one obtains a trascenden- 
tal equation for p* : 



2ae 



Bi r - 1 



+VI—VQ = 
(54) 

that one solves numerically for p* . This information pro- 
vides the necessary information to construct the solution 
in the whole domain as described below. The different 
regions are illustrated in Fig. |3] and the exact solution 
region by region is as follows. 

1. Region 1 is defined hy x — xq < tVshock, where the 
velocity of the shock is given by (|42| because the 
shock is traveling to the left: 



Vs ^vi - ai 



(r + i)P3 



1 



2pir 2r ' 

and the exact solution here reads 

Pexact Pi 1 

'^exact — 
Pexact Pi- 

2. There is no region 2. 

3. Region 3 is defined by the condition tVg < x — 
xo < tVcontact- Vcontact IS the characteristic value 
^ evaluated at this region: Vcontact = ^ V4 = 
V* . Using (p9| explicitly for the density and (|52)) 
for the velocity, the solution in this region reads 



Pexact 
'^exact 

Pexact 



P3, 
V3, 



Pi 



Pi(r-i) + P3(r + i) 
P3(r-i) + pi(r + i)- 



4. Region 4 is defined by the condition tVcontact < 
X — Xq < tVt , where the velocity of the tail of the 
rarefaction wave Vt is the third eigenvalue Q eval- 
uated at the region behind the tail Vt = W4 + 04. 

One uses (|53l) to calculate W4 and ([T3l) implies 
Pi/Pe — (Pi/pe)^ for a- constant value of K, which 
implies an expression for p^. The resulting exact 
solution is 



Pexact 
'^exact 

Pexact 



Pi, 
P6 



1/r 



Region 5 is a fan region defined by the condition 
tVt < X — Xq < tVh where the velocity of the head 
of the wave is again the third eigenvalue, but this 
time evaluated at the head Vh = vq + a^. One 
uses the expressions for a fan region of a rarefaction 
wave moving to the right (|27l28l29p to calculate the 
exact solution 



Pexact 



'^exact 



P6 



r - 1 



r + i a6(r + i) 



V6 - 



X — Xq 
1 



r + 1 

Pexact — P6 



-«6 + -^[T - 1)V6 H — 



F- 1 



F + 1 a6(F-|-l) 



V6 



X — Xq 
t 



6. Region 6 is defined by the condition tVh < x — Xq. 
The exact solution is given by the initial states at 
the right chamber. 



Pexact 
'^exact 
Pexact 



P6, 
P6- 



An example is shown in Fig. |3]for initial data in Table 
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3. Case 3: Rarefaction-Rarefaction 

A physical situation that provides this scenario is pL ~ 
PR, PL ~ PR and —vl = > 0. In this case both, 

regions 2 and 5 correspond to rarefaction waves. In this 
particular case since one of the rarefaction waves moves 
to the left and the other one to the right, we distinguish 
them using the labels for each of their parts. 

Again the contact discontinuity defines a relationship 
between velocity and pressure. In the present case, there 
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FIG. 3: Description of the relevant regions for the Shock- 
Rarefaction case. 



Again, once p* is calculated numerically, the solution in 
all the regions of the domain can be calculated as follows. 
The first implication is that = Pi = P* , and thus 
V3 and W4 can be calculated using (|55p and (|56|) . The 
different regions are illustrated in Fig. [5] 

1. Region 1 is defined by the condition x — xo< tVh,2, 
where Vh.2 is the velocity of the head of the wave 
moving to the left, and is obtained from the char- 
acteristic value of such rarefaction wave evaluated 
at the left interface, that is Vh.2 ~ vi ~ ai. In this 
region the gas has not affected the initial state on 
the left, then the solution is 



Pexact — Pi 1 
Vexact = Vi, 
Pexact — Pi- 





FIG. 4: Exact solution for the Shock-Rarefaction case at time 
t — 0.25 for the parameters in Table [T] 



is an expression for in terms of f 1 for a rarefaction 
wave moving to the left given by (jl9p and another one 
for V4 in terms of vq for a rarefaction wave moving to the 
right (1211): 



V3 = Vl 



2ai 

r - 1 

2a6 

r - 1 



- 1 



P6 



(55) 
(56) 



The condition = = v* at the contact discontinuity 
implies a trascendental equation for p* = = pa'- 



2a6 
r - 1 



1- - 



2ai 



- -1 



2. Region 2 is a fan region defined by the condition 
tVh,2 < X — xo < tVt^2, where the velocity of the 
tail Vt^2 is that of the state left behind by the wave, 
that is Vt^2 ^ V3- a^. 

The exact solution is that of a fan region of a rar- 
efaction wave moving to the left (j24l25l26p 



Pexact 



'^exact 



Pexact 



Pi 



2 r - 1 

+ 



r + i ai(r + i) 



Vl 



X — Xo 
t 



r + 1 
pi 



ai + ^{^- l)^'i + — - — 



r - 1 



r + i ai(r + i) 



Vl 



X — Xq 
t 



3. Region 3 is defined by the condition tVt^2 < x—xq < 
tVcontact- The velocity of the contact discontinuity 
is Vcontact = W3 = ^4 = according to the eigen- 
value ([3]). In this region p3 = p* and vj, ~ v* are 
already known from p* . Finally, the density is ob- 
tained from (fT3)l for an isentropic process like the 
rarefaction wave for a constant C on both sides of 
such wave as found in the previous two cases. Thus 
the solution is 



Pexact 
Vexact 

Pexact 



P3, 
V-i- 

Pi 



1/r 



+V1—VQ = 
(57) 



Region 4 is defined by the condition tVcontact < 
X — Xo < tVt^5, where the velocity of the tail of 
the wave moving to the right 14_5 is given by the 
eigenvalue (|H) evaluated at the state left behind 
the rarefaction wave moving to the right, that is 
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^,5 = ^'4 + 0.4, where again we point out that 114 = 
V* and P4 = p* are known once p* is calculated. 
The solution is obtained in the same way as for the 
previous region, but now the wave relates states in 
regions 4 and 6: 



Pexact 
'^exact 

Pexact 



Pi, 

Pe 



i/r 



5. Region 5 is defined by the condition tVt_5 < x — 
xo < Vh,5, where the velocity of the head of the 
wave moving to the right is Vh,5 = vq + uq, and the 
solution is obtained using the values of the state 
variables for the fan of a rarefaction wave moving 
to the right (|27I28I29|) : 




4 



FIG. 5: Description of the relevant regions for the 
Rarefaction-Rarefaction case. 



Pexact 



'^exact 



Pexact 



P6 



r - 1 



r + i a6(r + i) 



V6 



X — Xo 

t 



r + 1 

pe 



r- 1 



r + i a6(r + i) 



V6 



X — Xq 
t 



6. Finally, region 6 is defined by the condition Vh^5 < 
X — Xq. The exact solution is given by the initial 
values at the chamber on the right because in this 
region the gas has not been affected yet by the dy- 
namics of the gas: 



Pexact — Pe 7 
exact — ^61 
Pexact ~ Pe • 

An example is shown in Fig. [S]for initial data in Table 



4- Case 4-' Shock-Shock 

A physical situation that provides this scenario cor- 
responds to two streams colliding with opposite direc- 
tions. We choose in this case pl = pn, Pl = PR and 
—vl = +Vfi, < 0. In this case regions 2 and 5 are shock 
waves. 

Again the contact discontinuity defines a relationship 
between velocity and pressure. In the present case there 
is an expression for in terms of Vi for a shock-wave 
moving to the left given by (j4ip and another one for U4 
in terms of vq for a shock- wave moving to the right (1431) : 




FIG. 6: Exact solution for the Rarefaction-Rarefaction case 
at time t = 0.25 for the parameters in Tabled 



V3 = Vi- {p3 -pi)i 



Vi ^ ve + {p4 -pe)] 



P3 + Bi ' 



A. 



Pi + Be 



(58) 
(59) 



The condition W3 = W4 = v* at the contact discontinuity 
implies a trasccndcntal equation for p* = ps = Pi'. 



{P*-Pi)] 



Ai 



Ae 



p*+Bi 



-{P*-Pe)\l = 0. 

\ P ' + Bq 

Again, once p* is calculated numerically, the solution in 
all the regions of the domain can be calculated as follows. 
Immediately one has that pz ~ pa ^ p* and vt, and U4 
can be calculated using ((55)) and ([5^ . 
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In this particular case regions 2 and 5 reduce to lines. 
The solution in each region reads as follows and the re- 
gions are illustrated in Fig. [71 

1. Region 1 is defined by the condition x — xq< tT4,2, 
where the velocity of the shock moving to the 
left Vs^2 is given by and reads Vs^2 = wi — 

The solution there is that of 



Oi 



(r+i)P3 



2pir 



r-i 
2r ■ 



the initial values of the variables on the left cham- 
ber: 



2 






5 






/' ^ 






3 






1 






6 



Pexact Pi ■) 

'^exact — '^li 
Pexact ' Pi- 

2. There is no region 2. 

3. Region 3 is defined by the condition tVs^2 < x — 

Xq < tVcontact, where Vcontact = V3 = V4 = v"" . 

Once ([54)) is solved one can calculate all the re- 
quired information. Using ([58)) for ^3 and ([39)) for 
ps the solution in this region reads 



pexact 
^exact 

Pexact 



P3, 
V3 



Pi 



Pi(r-i)+p3(r + i) 
P3(r-i) + pi(r + i)- 



Region 4 is defined by tVcontact < x - xq < tVs^5, 
where the velocity of the shock moving to the 
right is given by (|45)) and reads Vs^^ ^ ve + 

Finally, using ^ 



„ (r+i)p4 



r-i 
2r ■ 



for and 

(^1)) for p4 the solution in this region reads 



FIG. 7: Description of the relevant regions for the Shock- 
Shock case. 



FIG. 8: Exact solution for the Shock-Shock case at time t = 
0.25 for the parameters in Table ID 



Pexact 
^exact 

Pexact 



P4, 
1)4, 



P6 



P6(r-i)+P4(r + i) 
P4(F-i)+p6(r + i)' 



5. There is no region 5. 

6. Finally region 6 is defined by the condition Vs.5 < 
X — xo- The exact solution is given by the initial 
values at the chamber on the right: 



Pexact — PG 1 
'^exact — ^65 
Pexact — PG- 



An example is shown in Fig. [5] 



III. RELATIVISTIC SHOCK TUBE 

First of all one needs to define a model for the gas. 
In our case we use the perfect fluid defined because it 
has no viscosity nor heat transfer, is shear free and is 
non-compressible. Such system is described by the stress 
energy tensor 



(61) 



where po is the rest mass density of a fluid element, 
its four velocity, p the pressure, h = \ + e + p/ po is 
the specific enthalpy and 77^'' are the components of the 
metric describing Minkowski space-time. 

The set of relativistic Euler equations is obtained from 
the local conservation of the rest mass and the local con- 
servation of the stress energy tensor of the fluid, which 
are respectively 
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(po«n,p = 0, 

where u'' = VK(l,u^,0,0) and = , ^ is the 

Lorentz factor and is the Eulerian velocity of the fluid 
elements. It is possible to arrange these equations as a 
flux balance set of equations as in the Newtonian case 



that all quantities describing the fluid depend on the vari- 
able ^ = [x — Xo)/t. In order to explore the change of 
all physical quantities along the straight line we define 
the useful change on the derivative operators 

dt = -]id^, d., = ^d^. (69) 

Using the advective derivative da = dt+vdx, we obtain 
the expressions 



dtvi + a,F(u) = 0, (62) 

where conservative variables are defined by u = 
{D,S^,t)'^ and the resulting fluxes are F = {Dv.Sv + 
p,S), where we assume that specifically v = v'^ and 
S ~ , since we are only considering one spatial di- 
mension. The conservative variables are defined in terms 
of the primitive ones as follows 

D = poW, 

T = = pohW-p. (63) 
The flux balance equations are explicitly: 



d^p - -DdaihWv), (70) 
dtP = Dda{hW), (71) 

where we have used the rest mass conservation law to 
simplify the expressions. From (|69p we obtain for the 
advective derivative da = — v)d/d^, for which we will 
use d := d/d^ from now on. With this in mind we obtain 
from (|70l7ip the differential equation 

{v - ClphW^dv + ^v)dp = 0. (72) 

On the other hand, the change of variable in (|64| from 
i, a; to ^ implies 



dtS^-dx{Sv + p) = Q, 
dtT + d^s = 0. 



(64) 
(65) 
(66) 



The eigenvalues of the Jacobian matrix of this system of 
equations are 



V ±Cs 
1 ± VCs 



(67) 



Each of the characteristic values ([67]) may correspond to 
eigenvectors with different properties exactly as in the 
Newtonian case, that is, A° corresponds to a contact dis- 
continuity, whereas the eigenvalues may correspond 
to rarefaction or shock waves. The shock tube problem 
in this case is defined as in the Newtonian case: 



u 



Ul, X < Xo 
Ur, X > Xo- 



(68) 



Next we describe the treatment of each of the wave or 
discontinuities that develop during the evolution. 



A. Rarefaction Waves 

Rarefaction waves are self-similar solutions of the fiow 
equations Q. They are self-similar solutions in the sense 



{v - ^)dp + pW^{l - v£,)dv = 0. 



(73) 



and from equations ([7^ and ([75)1 we obtain a relation 
between the density and pressure 



dp = h 



1 - < 



dp. 



(74) 



Since the process along ^ is isentropic [6[ the sound speed 
is = ^f^lsj which combined with the previous expres- 
sion implies the speed of sound 



l-< 



(75) 



Besides, we can find a useful expression for an isentropic 
process using p — Kp^ (we are using a politropic equa- 
tion of state). 



(76) 



From system (j67p we obtain the speed of sound in terms 
of the eigenvalues of the Jacobian matrix 



-{v- X+)/{l-vX+) ife = A+, 
{v- X-)/{l-vX-) if^^A-. 



(77) 
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Comparing with (1751) wc find that Cs{v, A"*") is the speed 
of sound for a rarefaction wave traveling to the right and 
Cs(u, A~) for a wave traveling to the left. 

According to this equation wc get from (|73)) that 



W^dv±—dp = 0. 



(78) 



Here the + sign refers to the wave traveling to the left 
and the — sign when it travels to the right. From this 
equation we obtain the Riemann invariant because this 
differential equation is valid along a straight line along 
the X — t plane, as long as it is not a shock. Integrating 
the first term of (|75)) we obatain 



where is 



Cs 



±2(r-i) 



-1/2 



(86) 



Equation (|85|) is valid only across straight lines aris- 
ing from the origin {xo,t = 0) and evolving along ^ = 
{x — xo)/t inside the rarefaction zone. For this family of 
straight lines the Riemann invariant is the same. This 
allows us to relate any two different states in the rarefac- 
tion zone, particularly we are going to take the states L 
and R as the states just next to the left and to the right 
from the rarefaction wave. 



1, 1 + v 
- in 

2 1-v 



± 



dp = constant. 



(79) 



In order to calculate the integral we use the definition 
of the sound speed and the polytropic equation of state 
p = Kp^ , from which we obtain 



clip) 



KTjT - l)p^-i 
r - 1 + ATpr-i 



(80) 



or in terms of the pressure instead of the density the 
speed of sound reads 



cUp) 



i~r ^ p_^ 
AT \K/ 



(81) 



1 



Conversely, if the speed of sound is known one can 
calculate the density using ([50]) : 



L 1 ^, R' 



Vl 



1 - 



(87) 



Assuming that when the wave is propagating to the 
left we account with information from the left state, we 
can calculate the velocity of the fiuid on the region at the 
right from the wave in terms of the state variables on the 
state at the left and 



VR 



{1 + vl)A+~{1~vl)A+ 
il + VL)A+ + {l~VL)A+- 



Analogously when the wave is moving to the right we 
expect to account with information on the state to the 
right. Then we can express the velocity on the left in 
terms of the variables on the state at the right and A~ 



{1+vr)Aj^ - {l-vii)A-^ 
(1 + vb)A^ + (1 - vb)A-^ 



(89) 



Then the integral can be written as 



(82) 



—dp = 
P 



1 



r - 1 



^dcs. (83) 
dcs 



Integrating by parts and using (|79p we find the useful 
constraint 



il 1 + "+ 1 

2 1 - u (F - 1)1/2 



In 



which in turn simplifies as follows 



1 

Y-v' 



-A^ 



constant^ 



constant^ 
(84) 

(85) 



1. The fan 

The fan is the region where the rarefaction takes place, 
propagating with velocity either A"*" if the wave is moving 
to the right or A~ when moving to the left. The fan will 
be bounded by two values of ^ corresponding to the head 
and the tail of the wave: 



VL.R ± c. 



(L.R) 



1 ± vc 



(L,R) ' 



Vr.L ± C, 



1 ± VC, 



(.R,L) ' 



(90) 
(91) 



where the — sign applies to waves traveling to the left and 
+ when the wave moves to the right. In order to con- 
struct the solution inside the fan, we use the constraint 
([57)1 . We have two cases according to the direction of the 
rarefaction wave. If the rarefaction wave travels to left 
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1 + ?^,+ 



l + VR 







(92) 



and solve the equation for vj^. When the rarefaction wave 
travels to right we use 



1 - Wfl 



0, 



(93) 



and solve the equation for v^. We calculate in each case 
using (j86p in the appropriate region 



A 



■^s,{L,R) 



^s,{L,R) 



±2(r-i)- 



(94) 



where the sound speed is given by ([751) ^^.nd ([77]) 



+ _ j_ V{L,R) - ^ 



(95) 



where the + sign is used when the wave moves to the left 
and — when moving to the right. Finally since we are in 
the rarefaction zone we can express a point (a;, t) with ^ ~ 
(x — xo)/t in (j95p and using this expression in (j94p and 
substituting into (|92|) or ()93p depending on the direction 
of propagation we finally obtain a trascendental equation 
for the velocity W( l,r)- We assume that if the wave moves 
to the left we know the variables on the state to the left L 
and ignore those of the state to the right R and viceversa. 
Then we look for a solution of when the wave moves 
to the left and of vr when moving to the right. Instead 
of looking for a closed solution to this equation we solve 
it numerically to obtain V(l,b.) assuming we know W(_r,l). 
Once W(L.i?,) is calculated we can substitute back, and 
using equation (|95|) obtain the sound speed; next, using 
([5^ obtain the density p; finally with the help of the EOS 
we can calculate the pressure p = Kp^ . This completes 
the solution in the fan region. 

The particular cases described later illustrate how to 
implement this procedure. 



B. Shock Waves 

Shocks rec^uirc the use of the rclativistic Rankinc- 
Hugoniot jump conditions [pQU^]n^ = and [T'^'^]n^ ~ 
across the shock 0], where n'^ ~ {—VsWs,Ws, 0,0) is 
a normal vector to the shock's front, Ws is the shock's 
Lorentz factor and Vg is the speed of the shock. Here 
we have used the notation [F] = Fl — Fr, where Fl 
and are the values of the function F at both sides 
of the shock's surface. These conditions reduce to the 
following system of equations, in terms of primitive and 
conservative variables, as 



Dlvl-Drvr = Vs{Dl~Dr), (96) 
Slvl +PL - (Srvji+pr) = Vs{Sl-Sr), (97) 
Sl~Sr = V,{tl-tr). (98) 

The subindices (L, R) represent two arbitrary states at 
left and at the right from the shock. These equations 
can be written in the reference rest frame of the shock 
by considering a Lorentz transformation, that is 



Dlvl = Drvr, 
Slvl+Pl = SrVr+pr, 
Sl = Sr, 



(99) 
(100) 
(101) 



where the hatted quantities are evaluated at the rest 
frame of the shock. Here V(^l.r) = T^Fl^^T^' ^(l.R) = 
P(l,r)Wlm., S^l.r) 



P{L,R.)h(^L,R)W^L,R.)^{L,R.) and 



Wi 



{L,R) 



From (j99p . we can introduce the invariant rclativistic 
mass flux across the shock as 



J = WsDl{Vs - Vl) = WsDr{Vs - vr), 



(102) 



where 



It is important to point out that 



when the shock moves to the right the mass flux is pos- 
itive j > 0, whereas when the shock moves to the left it 
has to be negative j < 0. 

Now, using the expression for the mass flux ()102|) into 
the Rankine-Hugoniot conditions (|96l [97l [98]) we can ob- 
tain the following system of equations in terms of a com- 
bination of primitive and conservative variables 



Vl - Vr 
PL - PR 
VlPL - VrPr 



i 




1 










Sl 


Sr. ' 




Dl 


Dr, 


^{ 


tl 


tr ' 


Ws V 


Dl 


Dr 



(103) 
(104) 
(105) 



Considering the shock is moving to the right and thus 
that the state R is known, we will write an expression for 
the velocity wi, in terms of the state variables R and also 
in terms of j, Vs and pl- In order to do this, we rewrite 
expressions (|104p and (|105p using the definitions for the 
conservative variables in terms of the primitive variables 
as follows 



Ws, , 
— 

Ws, . 

— {vlPl - VrPr) 

J 



hLWL-hRWR^, (106) 

VL 



HlWl - 
hRWR + 



PL 

PlWl 

PR 

PrWr' 



(107) 
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Subtracting these expressions and dividing by we get 



i 

HrWr 

PL 



VL 



vrPr 

PL 



1 



PR \ 

vlPl J 



(108) 



vr 

VL 



1 



PR 



1 



PlPrWr PlWl 



Inserting this into (|103p we finally obtain an expression 
for the velocity vl 



PLhLW^WliVs-VLf-pRhRW^Wl^iVs-VR)^ = -(pl-Pr) 

(114) 

As we can see from this equation, the definition of the 
conserved mass flux is present, then using equation (jl02p 
in this last equation, we obtain a useful expression for 
the square of the flux 



VL 



hRWRVR + ^{pL~ Pr) 



It-rWr + {pL - pr) 



j 



1 



PrWr 



(109) 



-{pl-Pr) 

(hL__ hR \ ' 

\PL PR J 



(115) 



When the shock moves to the left and the state L is 
known, the velocity on the state to the right is 



tiLWLVL + ^(pr-Pl) 



vr 



HlWl + [pR - PL 



)(^ 



PlWl 



(110) 



where the condition j < has to be satisfied. 

In order to obtain the shock velocity Vg , we start form 
the mass flux conservation across the shock (|102p , which 
relates the shock velocity with the mass flux. Substitut- 
ing Ws ~ l/-\/l — y/, it is possible to solve the resulting 
quadratic equation and obtain the two roots for the shock 
velocity 



V, 



p'rWI+P 



PlWlvL-Vf+fpj 

Plwl+f 



(111) 
(112) 



which correspond respectively to a shock moving to the 
right and to the left. The signs of the quadratic formula 
are chosen such that they are physically possible, that is, 
for the case of a shock moving to the right j > we use 
piip and for a shock moving to the left j < we use 

dmi) i. 

In order to solve completely the problem across the 
shock, we first express equation pOOp as 



where the positive root corresponds to a shock moving 
to the right whereas the negative root to a shock moving 
to the left. 



Another useful expression comes from equation (jlOip . 
which can be rewritten directly in the form 



hLWL = HrWr, 



which combined with equation (jllSp implies 



(116) 



hi 



hi 



[PL 



PR) 



hL 
PL 



hR 
PR 



(117) 



This last equation is commonly called the Tauh 's adiabat. 
Moreover equations (jllSp . ()116p and (|117p are known as 
relativistic Taub's junction conditions for shock waves 

m- 

Finally, in order to obtain the density pL and pressure 
PL for a shock moving to the right in terms of the vari- 
ables in the region to the right, we consider the definition 
of the specific internal enthalpy and that the fiuid obeys 
and ideal gas equation of state. With these assumptions 
equation (|117p can be rewritten in the form 



PLhL{Vs - VlY 



Vs'vl 



PRhRjVs - vrY 
1 - 1/2 - vl + 



s ""R 



- {PL - PR 

Considering that WsW(^l.r) 



(113) 



the last equation takes the follow- 



ing form 



— [pl{2<j-1)+ pr] + ^[pl{a-l)+ plPr] = 

PL Pl 

— [pR{2a - 1)+pl] + ■^[Pr{<J - I) + plPrKUS) 
PR Pr 



where a — y&j- The solution for the quadratic equation 
reads 
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1 

PL 
1 
PR 



-[pL{2a - 1)+PR\ ± ^[pl{2(J - l)+PR? + 4ag[pi(g - 1)+PlPr] 

2a[pl{(j-l)+pLPR] 
-[pBi^a - 1) + Pl] ± v/M2a - 1) + p^Y + 4Cflcr[p^(cr - 1) + p^pl] 



2a[pl{a-l) 



PrPl\ 



(119) 
(120) 



where Cl = ■j^[pR{2a-\) 



Pr] - 



-[pl{f^ - 



- 1) +PlPr] 

1) +PflPL] 



andCii = ^bL(2a-l) _ . , 
A physically acceptable solution requires p > 0, which 
restricts the sign to be positive one in both cases. 



C. Contact Wave 

The equations describing the jump conditions 
(|96I97I98P admit the solution using Vs = vr = vl = 
X° ~ Vcontact where vr and vl are the values of the ve- 
locity of the fluid at the right and at the left from the 
contact discontinuity. This represents the contact wave 
traveling along the line x — xq = A°t. 

Then ([M]) is trivial and reads 



_ {l + v,)A+ -{l^v,)A+ 
{l+vi)A+ + {l-vi)A+' 

where according to ([M)) 



A 



(1,3) 



s,(l,3) 



.(1,3) 



(122) 



(123) 



Herec+i := Cs{pi) = V'lW(PiM> = l+ ^Jr-i) 



and 



-8,3 



CsiPa) is given by equation 



(Sl - Sr)Vs +PL-PR = {Sl - Sr)Vs , 



(121) 



which implies pR = pl and equation (|98p is satisfied. 

We are now in the position of analyzing each of the 
possible combinations of shock and rarefaction waves in 
a Riemann problem. We then proceed in the same way 
as in the Newtonian case studying each combination. 



D. The four different cases 



fe) 



\ r-1 

\ KT 



K 



r-i /PS ^ 

KT \K) 



1 



Pl_ 

pV 



(124) 



where we remind the reader that in the rarefaction region 
the polytopic constant remains the same during the pro- 
cess, that is, it is the same in regions 1, 2 and 3. On the 
other hand the velocity of the gas in region 4 corresponds 
to the velocity on the state at the left of a shock moving 
to the right (fT09)) 



In what follows, as we did for the Newtonian case, we 
present the four combinations of rarefaction and shock 
waves associated to the relativistic Riemann problem. 
We illustrate each case with a particular set of parame- 
ters contained in Table HIl 

1. Case 1: Rarefaction-Shock 

The contact wave conditions are ~ V4 = v* and = 
P4 = p* ■ The velocity in region 3 is given by equation 
((88)) that provides the velocity on the state at the right 
from a rarefaction wave moving to the left: 



W6 + (P4-P6)(^ + ^)' 



where Ws,5 = l/yl — V^^ is the Lorcntz factor of the 

shock, where we use the subindex 5 in order to denote 
the shock occurring in region 5. In order to obtain V4 in 
terms of p4 we need to perform the following steps: 

• The rest mass density p4 is given in terms of p^ and 
other known information can be expressed using 
(dm) as 



J_ ^ -b4(2g - l)+P6] + V[Pi(2a - l)+pe]^ + ^Mpjia - 1) + p^pe] ^^^g) 
P4 2(t[pI{<t - 1) + p4pe] 

C4 = —[P6i2<J-l)+p4] + -^[pl{a-l)+p4P6], where cr = -^. (127) 

P6 Pi r-1 
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Case 


PL 


PR 


VL 


VR 


PL 


PR 


Raref act ion- Sho ck 


13.33 











10 


1 


Shock-Rarefaction 





13.33 


0.0 


0.0 


1 


10 


Rarefaction-Rarefaction 


0.05 


-0.05 


-0.2 


0.2 


0.1 


0.1 


Shock-Shock 


3.333e-9 


-3.333e-9 


0.999999 


0.999999 


0.001 


0.001 



TABLE II: Initial data for the four different cases. We choose the spatial domain to be s £ [0, 1] and the location of the 
membrane at xo = 0.5. In all cases we use F = 4/3. 



Once p4 is given in terms of p4 it is possible to 
compute the enthalpy in region A as = 1 + a^. 

Then equation (|115p reads 



j2 ^ (Pi - Pe) 



h4_ 

Pi 



P6 



(128) 



where ha ~ I + a—. Something to remember here 
is the fact that as the shock moves to the right, we 
consider j to be the positive square root. 

Once j is obtained, the shock velocity can be found 
from expression (|llip as 



s,5 



plWive + \j\Vf 



Pi 



f + Piwi 

• Finally one calculates Ws,5 = 



(129) 



and in this 



way W4 in terms of p4 and the known state in region 
6 using (fT25|) . 

According to the contact discontinuity condition = 
V4 = V* , we equate (|122p and (|125p and obtain a tran- 
scendental equation for p* : 



(l + vi)A+ + (l-wi)A+(p*) 



heWe + ip* -pe) 



peWe 



= 0, (130) 



which has to be solved using a root finder. 

Once this equation is solved, and p4 are automati- 
cally known and ^3 and V4 can be calculated using (|122p 
and (|125p , respectively. It is possible to calculate pa using 
the fact that in the rarefaction zone the process is adi- 
abatic and then p3 — Pi{pz/piY^^ ■ On the other hand 
we can also calculate using (|126l) . With this informa- 
tion it is already possible to construct the solution in the 
whole domain. 

Up to this point we account with the known initial 
states {pi,vi,pi) and {pq,vq, po), the solution in regions 
3 and 4 given by (pa, 1^3,^3) and (p4,W4,P4), and 14,5 
which represents the velocity of propagation of the shock 
5. The exact solution region by region is described next. 



Region 1 is defined by the condition x — xq < t^h, 
where according to (PU)) is the velocity of the 
head of the rarefaction wave traveling to the left 



The values of the physical variables 



l — ViCs^ 

are known from the initial conditions: 



Pexact — Pl-} 
'^exact — ^^1 1 
Pexact — Pl- 



(131) 
(132) 
(133) 



2. Region 2 is defined by the condition t£^h < x — xq < 
tS^t, where according to ((9T|) is the characteristic 
value again, but this time evaluated at the tail of 
the rarefaction wave, that is £t = "^"'^'■^ . In order 

to compute V2 we use ([5^ 



\±^At-\±^Ativ2)^0 (134) 
considering equations ([7S)) . (|M| and (^5)) as follows 



1 +2{T-l)-^'^ 



^(1,2) 




1 - V2^ 



where ^ = (x — xo)/t. In this way, equation (|134p 
is transcendental and has to be solved equivalcntly 
for V2 or for cj'j using a root finder for each point of 
region 2. We recommend solving for cj^j and then 
construct V2 using (|137p . Finally we calculate p2 
using equation (|82|): 



P2 



KT 



1 

r-1 



Pl 



(138) 
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Finally we obtain p2 using the fact that in the pro- 
cess K is constant 




(139) 



3. Region 3 is defined by the condition t^t < x ~ xq < 
tVcontact, where Vcontact = Aq = U3 = V4- The 
solution there reads 

Pexact = P3, (140) 
Vexact = V3, (141) 
Pexact = P3- (142) 

4. Region 4 is defined by the condition tVcontact < 
X — xq < tVs^5, where Vs^5 is given by (|129|) and 
explicitly 



2. Case 2: Shock-Rarefaction 

This is pretty much the previous case, except that one 
has to be careful at using the correct signs and conditions. 
We then start again with the contact wave conditions 
V3 = Vi — V* and P3 = Pi = p* ■ The velocity of the gas 
in region 3 corresponds to the velocity on the state at the 
right from a shock moving to the left (|110[) 



/^iTVii^i + ^(P3-Pi) 

h,w,+{p,-p,)[^+^y 



where Ws.2 ~ Vy 1 ^ ^s2 the Lorentz factor of the 

shock. In order to obtain in terms of py, and other 

known information we need to perform the following 
steps: 



Pexact — Pa, 
Vexact — ^4: 
Pexact = Pi- 



(143) 
(144) 
(145) 



5. There is no region 5. Only the shock traveling with 
speed 14,5- 

6. Region 6 is defined by tVs^^ < x — xq. In this region 
the solution is simply 



Pexact — P&^ 
exact — ^6; 
pexact P6 • 



(146) 
(147) 
(148) 



As an example we show in Fig. |9] the primitive vari- 
ables dX t = 0.35, for the initial parameters in Table HIl 




FIG. 9: Exact solution for the Rarefaction-Shock case at time 
t — 0.35 for the parameters in Table HIl 



The rest mass density is given in terms of p^ using 
the expression (jl20p as 



J_ ^ -[P3{2a - l)+Pi] + V[P3{2(J - 1)+Pl? + 4C3g[P3(^ - 1) +P3Pl] 

P3 2a[plia - 1) + p3Pi] ' ^ ' 

Cs = — bi(2cr - 1) +P3] + -^[pi(cr - 1) +p3pi], where cr = -. (151) 

Pi Pi r - 1 



• Once p3 is given in terms of p^ it is possible to 
compute the enthalpy in region 3 as ^,3 = 1 + tr^. 



• Then from equation pi5|) we obtain 



■2^ {P3-Pi) 



h3 

P3 



hi 
Pi 



(152) 



where hi = 1 + . As the shock is moving to 
^ pi ° 
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the left we consider the negative root of the above 
expression for j. 

• Once j is obtained, the shock velocity can be found 
from expression ()112p in terms of as 



V, 



s,2 



• Finally one calculates Ws^2 = 



(153) 



and in this 



way W3 in terms of and the known state in region 
1 using (|149p . 

The velocity in region 4 is given by equation (|89|) that 
provides the velocity on the state at the left from a rar- 
efaction wave moving to the right: 



1. Region 1 is defined by the condition x — x^ < tVs.2j 
where Vs. 2 is given by (|153p and the solution there 
is that of the initial state on the left chamber 



Pexact 
'^exact 
Pexact 



Pi, 
Vl, 
Pi- 



(158) 
(159) 
(160) 



2. There is no region 2. Only the shock traveling with 
speed Vs,2- 

3. Region 3 is defined by the condition tVs,2 < x — 
xq < tVcontact, where Vcontact = Aq = W3 = W4- The 
solution is 



1^4 - — — — 

(1 + (76)^6 +(1-1^6)^4 

where following ([M)) 



(154) 



A 



,(4,6) 



,(4,6) 



-2(r-l) 



Here ^ 



(4,6) 



Cs{pi) is given by equation 



(155) 



CsiPe) = \/^P6/{P6he), he = l + 



PeT 

P6(r-i) 



and 



\ 



r- 1 



K 



r-i 



P4, ' 



El 
pV 



(156) 



because K is the same in regions 4 and 6. 

We obtain a transcendental equation for p* using the 
contact discontinuity condition ^3 = 774 = u* , and equate 
([T49| and ((T54l) : 



Pexact — PSi 
exact — ^3: 
Pexact — PZ- 



(161) 
(162) 
(163) 



4. Region 4 is defined by the condition tVcontact < 
x — xo < t^t, where according to is the 

characteristic value again, but this time evaluated 
at the tail of the rarefaction wave, that is = 



1 + ViCs.i 



The solution in this region is 



Pexact — Pa, 
Vexact — "^i, 
Pexact = Pa- 



(164) 
(165) 
(166) 



5. Region 5 is defined by the condition tS^t < x — xq < 
tS^h, where according to ((90|) ^/^ is the velocity of the 
head of the rarefaction wave traveling to the right 



V6 + Cs,6 



In order to compute W5 we use ([M)) 



{l + ve)A^ -{l-ve)A^{p*) 
il + ve)Ae +{l-ve)A^{p*) 
hiWivi + ^^{p* -Pi) 

hiW, + {p* - p,) + ^) 



0, (157) 



which has to be solved using a root finder. 

Once this equation is solved, ps and p4 are automati- 
cally known and f 3 and V4 can be calculated using (|149p 
and (jl54p . respectively. It is possible to calculate p4 us- 
ing the fact that in the rarefaction zone the process is 
adiabatic and then = peiPi/PeV^^ - We can also cal- 
culate p3 using (|150p . With this information it is already 
possible to construct the solution in the whole domain. 

Up to this point we have the known initial states 
{pi,vi,pi) and (j)6,ve, pe), the solution in regions 3 and 
4 given by (p3, t'3, ps) and (p4, 174, P4), and K,2 which rep- 
resents the velocity of propagation of the shock 2. The 
exact solution region by region is described next. 



-Ar^ - z Ar, [V5) = 0, 



1 - ve 



1- V5 



(167) 



whoch requires the information in (j76p . (|94p and 

mi: 



A5.6) 



^s,6 



,(5,6) 



,(5,6) J 



-2(r-i)-i/2 



, , , P6 



y Pehg 

V5 - £. 

l-f5^ 



V5 = 



P6 vr-1 



, (168) 

(169) 
(170) 



where = {x — Xo)/t. In this way, equation (|167p 
is transcendental and has to be solved equivalently 
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for ^5 or for cjg using a root finder for each point of 
region 5. We recommend solving for cj^ and then 
construct using (|170p . Finally we calculate ps 
using equation (|82|) : 



P5 



AT 



(171) 



since K is the same in regions 5 and 6, and by the 
same reason we obtain ps using 



P5 =P6 



(172) 



6. Region 6 is defined by t^/^ < x — xq. In this region 
the solution is simply 



Pexact — PGi 
exact — ^6: 
Pexact — PQ ■ 



(173) 
(174) 
(175) 



As an example we show in Fig. [TU]the primitive vari- 
ables aX t = 0.35 for the initial data in Table [TTl 



^ {l + vi)A+ -{l-vi)A+ 
{l+vMt + {l-vi)At' 

where according to ((94|) 



(176) 



A 



(1,3) 



,(1,3) 



-s,(l,3)J 



(177) 



Here 



-s,3 



c^i Cs(pi) = y/rpi/{pihi), hi = l + ^^|p^ and 



CsiPa) is given by equation (|8T 



^3(^3) = 



\ 



F- 1 



r-i /ps ^ 



1 



K = ^. (178) 

Pi 



On the other hand the velocity of the gas in region 4 
corresponds to the velocity on the state at the left of a 
rarefaction wave moving to the right 



(179) 



(1 + 1-6)^6 -(i- i'6)^r 

V4 = — 

(1 + 1-6)^16 +(l-t'6)^4 

where according to ([M)) 




FIG. 10: Exact solution for the Shock-Rarefaction case at 
time t = 0.35 for the parameters in Table ITIl 



3. Case 3: Rarefaction-Rarefaction 

In this case the transcendental equation for the pres- 
sure at the contact discontinuity is given again by the 
condition ^3 = where both velocities are constructed 
using the information of the unknown state aside rarefac- 
tion waves. The velocity in region 3 is given by equation 
([55)1 for the velocity on the state at the right from a rar- 
efaction wave moving to the left: 



A 



(4,6) 



"s,(4,6) 



-s,(4,6) 



-2(r-i)-i/2 



(180) 



and the speed of sound in region 4 is given by 



1(^4) = 



F- 1 



1 



K = ^. (181) 



Then using the contact discontinuity condition U3 — 
V4 — V*, we equate (|176p and (|179p and obtain a tran- 
scendental equation for p* : 

{l + vi)At-{l-vMt{p*) 



il + vi)A+ + {l-vi)A+{p*) 
{l + ve)A^ -il-ve)A^{p*) 
{l + ve)A^ + {l-ve)A^{p*) 



(182) 



which has to be solved using a root finder. 

Once this equation is solved, p3 and p^ are automati- 
cally known and ^3 and V4 can be calculated using (|176p 
and (|179p . respectively. As in the previous two cases, it 
is possible to calculate ps and p4 using the fact that in 
the rarefaction zone the process is adiabatic and then 
P3 = PiiPa/pi)^^^ and p4 = PQ{p4/peY^^ ■ Thus we 
have the known initial states (pi, pi), {pq,vq, po) and 



21 



the solution in regions 3 and 4 given by (pail's, pa) and 
(p4, P4). The solution in each of the fan regions aside 
the rarefaction zones has to be constructed in terms of 
the position and time = {x — xq) /t as described below 
for regions 2 and 5. 

1. Region one is defined by the condition x—xq < t^h2, 
where according to (|90p £^^2 is the velocity of the 
head of the rarefaction wave traveling to the left 
^h2 = il^^c '^i ■ values of the physical variables 
are known from the initial conditions; 



Pexact — Pi 1 
^exact — ^1 1 
Pexact Pi- 



(183) 
(184) 
(185) 



Region 2 is defined by the condition t£h2 < x^xq < 
t£t2, where according to (|9T|) ^42 is the characteristic 
value again, but this time evaluated at the tail of 
the rarefaction wave, that is £t2 = . In 

order to compute V2 we use (j92| 



3. Region 3 is defined by the condition t£t2 < x — xq < 

tVcontact, where Vcontact = Xo = V3 = Vi- The 

solution there reads 



Pexact — PS , 
^exact — ^3 , 
Pexact — P3- 



(192) 
(193) 
(194) 



Region 4 is defined by the condition tVcontact < 
X — xq < i^ts, where ^ts is the third characteristic 
value calculated at the tail of rarefaction moving to 
the right, and according to (|9T|) £,t5 = i'+t>ic \ ■ 
this region thus 



Pexact = P4, 
Vexact — ^^4, 
Pexact = Pa- 



(195) 
(196) 
(197) 



5. Region 5 is defined by the condition tS^tb < x — xq < 



t^hb- where ^,,5 



l+vac 

order to compute we use ([M)) 



according to ((90|) . In 



1 + 



1 + 1^2 



I -Vi 1 - V2 

where using dTB]), dM]) and ([^5]) 



A+(z;2)=0, 



(186) 



0, 



1 - W5 1-^6 

where according to (|76)) . ([94]) and ([95 



(198) 



4+ 

^(1,2) 



1 +2(r-i 



■(1,2) 



-s,(l,2)J 



,-1/2 



1-^2? 



. hi = l 



Pi 
Pi 



«2 



r - 1 



^s,2 



1 + <2^ 



(187) 

(188) 
(189) 



where £ = {x — Xo)/t. In this way, equation (|186p 
is transcendental and has to be solved equivalently 
for V2 or for using a root finder for each point 

of region 2. We solve for and construct V2 using 
(|189p . Finally we calculate p2 using equation ([82]) : 



i ^- Pi 

P2 = ^, A = — . 

■ - — Pi 



1 

r-1 



(190) 



Finally we obtain p2 using 



P2 =Pi 



(191) 



"(5,6) 



1 -2(r-i)-i/^ 



,(5,6) 



■,(5,6) J 



/ rp6 

Pehe 

«5 - C 
1 - ^'5^ 



P6 



/J6 = 1 + - 



P6 vr- 1 

" 1 - c-^e 



(199) 

(200) 
(201) 



where £, = {x — XQ)/t. Again (|198p is a transcen- 
dental equation either for U5 or for cjg. Once cjg 
has been calculated use (|20ip to construct U5 or di- 
rectly solve (jl98p for 1)5. It is possible to calculate 
P5 using ([82]): 



1 



P5 



and finally the pressure 

P5 =P6 



, (202) 



(203) 
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6. Region 6 is defined by t(^h5 < x ^ 
the solution is simply 



a;o. In this region 



4- Shock-Shock 



Pexact — P6 1 
f^exact — ^6: 
Pexact — PG- 



(204) 
(205) 
(206) 



We proceed as always, by establishing a relationship 
between the velocity in regions 3 and 4. We start by 
expressing V3 as the velocity of the gas on a region at the 
right from a shock moving to the left, that is, according 

to (inni) 



As an example we show in Fig. [TT]tlie primitive vari- 
ables at t = 0.25, for the initial parameters in Table llll 



hiWiVi 



V3 



Ws. 

J2 



■iP3-Pl) 



hiWi + (p3 -pi) 



32 



PiWi 



(207) 



V^2 is the Lorcntz factor of the 



where Ws^2 ~ V 
shock moving to the left. In this particular case we dis- 
tinguish between the two values of j depending using the 
subindices 2 and 5. In order to obtain 1)3 in terms of 
we can proceed following these steps: 



FIG. 11: Exact solution for the Rarefaction- Rarefaction case 
at time t = 0.25 for the parameters in Table [III 



• The rest mass density is given in terms of pa using 
the expression (|120p as 



]_ ^ -b3(2g - l)+pi] + ^[p3{2a-l)+pi]^ +AC3<j[pI{<J-1)+P3Pi] ,208) 
P3 2a[pl{a - 1) + p3Pi] 

Cs = —[pi{2a-l)+p3] + ^[pl{<J-l)+p3Pi], where f7 = -^. (209) 

Pi Pi r - 1 



Once p3 is given in terms of P3 it is possible to 
compute enthalpy in region 3 as ^,3 = 1 -I- cr 



P3 



Then from equation ()115|) we obtain 



• Finally we calculate Ws 



and thus V3 in 



terms of p3 and the known state in region 1 using 

HMD. 



(P3-Pl) 



h3 

P3 



hi 
Pi 



Using the information of the shock moving to the right 
(210) obtain the velocity at the left from the shock, that is 

V4 using (|109p 



where we choose j2 to be the negative root since 
the shock is moving to the left; here hi = 1 + a 



pi 



• Once j2 is obtained, the shock velocity can be found 
from expression pi2p in terms of p^ as 



Vs. 



pIW?vi-\j2\VjIT7i 



j1 



pIWi' 



(211) 



V4 



heWGVe + ^iP4 



Pa) 



heWe + (p4 ~Pe 



)( 



35 



sWe J 



(212) 



where Ws.5 = 1/ — V^r-, is the Lorcntz factor of the 

shock. In order to obtain 1)4 in terms of p4 wc need to 
perform the following steps: 
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• The rest mass density is given in terms of using 



the expression (|119l) as 



J 



1 

C4 



-[Pi{2a - 1) +Pg] + ^[p4(2ct - 1)+P6? + 4C4fTbl(tT - 1) + p^Pe] 



2a[pl{a - 1) + pipa] 



1 

P6 



[pe{2a - I) + Pi] + —[plia - 1) 



r 



-P4P6\, 



where cr — 



r- 1 



(213) 
(214) 



• Once p4 is given in terms of p4, we are able to 
compute enthalpy in region 4 as 



1 + 

Pi 



• Then equation (|115p reads 



J5 



h4 
Pi 



he 
P6 



(215) 



4 given by (^3,^3,^3) and {p4,V4,P4), together with Vs,2 
and Vs^5 which represent the velocities of propagation of 
the shocks. 

1. Region 1 is defined by the condition x — xq< tVs^2, 
where the velocity of the shock is (|21ip . The solu- 
tion there is that of the initial values of the variables 
on the left chamber: 



here he = I + a—. In this case, since the shock 
is moving to the right we choose the js to be the 
positive root. 

• Once j5 is obtained, the shock velocity can be found 
from expression (jllip 



s,5 



■pI 



j1 



piwi 



(216) 



and in this 



• Finally we calculate Ws 5 = , ^ 

way we can obtain V4 in terms of p^ with (j212p and 
the known state in region 6. 

According to the contact discontinuity condition ^3 = 
V4 = V* , we equate (|207p and (|212p and obtain a tran- 
scendental equation for p* : 



Pexact — Pli 
'^exact — ^^1 1 
Pexact — Pi- 

2. There is no region 2, only the shock wave traveling 
at speed Vs.2- 

3. Region 3 is defined by the condition tl4,2 < x ~ 
xo < tVcontact, where the velocity of the contact 
discontinuity is the characteristic value A° = u eval- 
uated in this region Vcontact ^ V3 ^ V4 ^ v* . 



Pexact 
'^exact 
Pexact 



P3, 
P3- 



4. Region 4 is defined by tVcontact < x 
and the solution is 



2^0 < tVs.B 



hiWivi + ^{p* -Pi) 
h,W, + {p*-p,)\^ + j^) 



heWe 



+ iP*-P^){^ + Jv;) 



0, (217) 



which has to be solved using a root finder. 

Once this equation is solved, p^ and p^ are automati- 
cally known, and V3 and V4 can be calculated using ([207]) 
and (|212p . respectively. It is possible to calculate ps and 
P4 using (|208p and (|213p . respectively. With this infor- 
mation it is already possible to construct the solution in 
the whole domain. 

Up to this point we have the known initial states 
(pi,wi,pi) and {pQ,ve, pq), the solution in regions 3 and 



Pexact 
'^exact 
Pexact 



P4, 
V4, 

P4- 



5. There is no region 5, only the shock wave traveling 
with speed Vs,^. 

6. Finally region 6 is defined by the condition Vs.5 < 
X — xq. The exact solution is given by the initial 
values at the chamber at the right: 



Pexact — PGi 
^ exact — ^6: 
Pexact — PQ- 
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As an example we show in Fig. [12] the primitive vari- 
ables a.t t = 0.55, for the initial parameters in Table HIl 





in the newtonian and relativistic regimes, which accord- 
ing to our experience is not presented in a straightforward 
enough recipe in literature. 

The contents in this article can be used in various 
manners, specially to: i) test numerical solutions of the 
Newtonian Riemann problem in basic courses of hydro- 
dynamics, ii) test numerical implementations of codes 
solving hydrodynamical relativistic equations, iii) under- 
stand the different properties of the propagation of the 
different type of waves developing in a gas and the differ- 
ent conditions on the hydrodynamical variables in each 
case. 

It is also helpful because with our approach it is pos- 
sible to straightforwardly implement the exact solution, 
and this will save some time to a student starting a career 
in astrophysics involving hydrodynamical processes. 



FIG. 12; Exact solution for a Shock-Shock case at time t 
for the parameters in Table HIl 
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